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Abstract 

We  develop  our  previous  work  on  hyperdistributions  into  a  formally  well-defined 
transform  which  may  be  applied  to  images,  the  hyperdistribution  transform  (IIDT). 
The  HDT  has  many  properties  in  common  with  conventional  orthogonal  transforms  of 
signals,  such  as  the  Fast  Fourier  Transform,  which  suggests  the  possiblity  of  develop¬ 
ing  a  fast  algorithm  for  the  IIDT.  Presently,  we  have  formulated  the  HDT  in  matrix 
language,  which  permits  a  reasonably  efficient  computational  approach  to  calculating 
the  HDT  of  an  image.  We  then  apply  the  HDT  to  image  compression  by  representing 
the  image  as  a  truncated  HDT  expansion  and  reconstructing  the  image  from  the  trun¬ 
cated  HDT  expansion.  The  compression  ratio  is  measured  in  terms  of  the  number  of 
bits  in  the  truncated  HDT  expansion  compared  to  the  number  of  bits  in  the  original 
image.  Test  cases  involving  both  synthetic  and  natural  images  are  considered.  Good 
quality  reconstructions  of  natural  images  are  obtained  with  compression  ratios  of  4:1 
and  recognizable  images  are  obtained  with  compression  ratios  of  16:1.  It  was  not  nec¬ 
essary  to  segment  the  images  into  sub-images.  Substantial  further  improvements  in  the 
performance  of  HDT  compression  may  be  obtained  by  employing  image  segmentation 
and  other  standard  techniques  for  transform-based  image  compression  algorithms. 
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1  Introduction 


Even  modest-sized  images  involve  large  amounts  of  data.  A  512  x  512  image  in  which  each 
pixel  may  take  on  256  gray  scale  values  will  require  1/4  Mbyte  of  computer  memory  for 
storage.  (1  Mbyte  =  1,048,576  =  bytes)  Accordingly,  storage  of  large  numbers  of  images 
places  large  demands  on  computer  memories.  Transmission  of  images  puts  correspondingly 
large  demands  on  channel  bandwidth  and/or  time  required  for  transmission  of  images.  For 
these  reasons,  among  others,  ther  are  considerable  economic  and  practical  pressures  to  de¬ 
velop  representations  of  images  which  require  smaller  numbers  of  bits. 

Fundamentally,  compression  algorithms  exploit  the  redundancy  (in  an  information  theo¬ 
retic  sense)  of  the  pixel  intensities  in  an  image.  This  redundancy  is  reflected,  for  example, 
in  the  statistical  dependence  of  these  pixel  values.  A  compressed  image  will  require  fewer 
bits  than  the  original  image.  Generally  speaking,  image  compression  techniques  fall  into 
two  categories^,  those  which  depend  on  coding  techniques,  such  as  Huffman  coding,  and 
those  which  depend  on  transform  techniques,  such  as  the  discrete  cosine  transform  (DCT) 
and  the  hyperdistribution-based  methods  which  we  have  developed  as  the  subject  of  this 
proposed  research  program.  We  shall  concentrate  attention  on  the  class  of  transform-based 
compression  algorithms. 

Exact  reconstruction  of  an  image  may  be  obtained  by  lossless  algorithms,  usually  at 
some  cost  in  the  compression  ratio  achievable  (typically  ~  10).  Applications  which  can 
tolerate  approximate  reconstructions  of  images  may  utilize  lossy  algorithms  for  compres¬ 
sion/decompression.  These  lossy  algorithms  may  be  particularly  appropriate  when  the  least 
significant  bits  (LSB’s)  of  an  image  are  corrupted  by  noise  and  thus  do  not  represent  de¬ 
sirable  data.  Lossy  compression/decompression  algorithms  can  achieve  higher  compression 
ratios  (40  or  more). 

A  number  of  techniques  have  been  developed  for  image  compression/decompression,  such 
as  the  DCT  and  techniques  based  on  affine  transformations  and  fractals.  Such  techniques  are 
computationally  intensive  and,  particularly  for  the  case  of  the  DCT,  may  require  implemen¬ 
tation  on  special  hardware  to  achieve  compression/decompression  of  images  in  reasonable 
times.  We  note  that  the  marketplace  for  commercial  applications  in  image  compression  is, 
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in  fact,  moving  in  the  direction  of  special  purpose  chips  for  rapid  execution  of  particular 
algorithms. 

A  compression/decompression  algorithm  which  may  be  used  in  conjunction  with  other 
image  processing  algorithms,  such  as  spatial  fdtering,  edge-enhancement,  or  deconvolution 
techniques  may  have  special  advantages.  In  our  Phase  I  research  we  have  developed  such  an 
approach  based  on  newly  developed  mathematical  called  hyperdistributions.  Hyperdistribu¬ 
tions  are  theoretically  attractive  because  they  are  an  algebraic  field  in  which  convolutions 
are  the  multiplicative  operation.  Thus,  problems  for  which  the  calculation  of  the  convolu¬ 
tion  inverse  (deconvolution)  is  not  well-posed  by  conventional  Fourier  transform  techniques, 
may  be  solved  uniquely  with  hyperdistribution  techniques.  This  property  also  is  reflected 
in  the  comparative  computational  stability  and  efRciency  of  deconvolution  computations  by 
hyperdistribution  techniques. 

On  the  basis  of  connections  with  the  problem  of  moments  and  formal  analogies  with 
moment  expansions,  we  believed  that  hyperdistribution  expansion  techniques  would  have  a 
significant  utility  for  problems  in  image  compression.  In  practice,  this  would  mean  that  an 
image  would  be  represented  in  a  two-dimensional  hyperdistribution  expansion  which  would 
be  truncated  at  a  finite  number  of  terms.  The  image  could  then  be  reconstructed  to  some 
level  of  precision  from  its  hyperdistribution  expansion,  in  a  fashion  analogous  to  that  we  have 
demonstrated  for  using  hyperdistributions  for  computing  global  approximations  of  functions. 
We  expected  that  for  many  images,  the  number  of  bits  required  for  the  hyperdistribution 
expansion  would  be  significantly  smaller  than  the  number  of  bits  in  the  original  image,  i.e. 
that  the  hyperdistribution  representation  of  an  image  is  a  valid  compressed  representation 
of  the  image. 

Our  Phase  I  research  effort  developed  the  mathematical  structure  required  for  effectively 
carrying  out  these  algorithms,  and  then  conclusively  demonstrated  the  validity  of  this  point 
of  view.  We  include  in  this  proposal  some  of  the  first  demonstrations  of  image  compres¬ 
sion/decompression  using  hyperdistributions.  We  are  able  to  achieve  compression  ratios  of 
16:1  routinely  even  with  algorithms  at  the  present  crude  stage  of  development.  We  believe 
that  by  fine-tuning  these  algorithms  and  combining  them  with  other  conventional  image  pro¬ 
cessing  algorithms,  that  higher  compression  ratios  and  improved  fidelity  of  the  reconstructed 
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image  is  possible.  Our  proposed  Phase  II  effort  will  lead  to  hyperdistribution  algorithms  for 
image  compression/decompression  competitive  with  other  entries  in  the  commercial  arena. 

We  establish  an  explicit  method  for  image  compression  and  reconstruction  using  hyper¬ 
distribution  theory.  The  point  function  approximation  for  the  hyperdistribution  expansion, 
the  Rodriguez  hyperdistribution  expansion,  is  used  to  formulate  a  hyperdistribution  trans¬ 
form  (HDT).  The  transform  includes  an  adjustable  parameter  which  is  used  to  vary  the 
shapes  of  the  reconstructing  wave  and  optimize  the  reconstruction  performance.  A  matrix 
representation  of  the  transform  similar  to  other  methods  is  derived.  This  establishes  algo¬ 
rithms  for  carrying  out  HDT’s  analogous  to  those  for  conventional  orthogonal  transforms 
such  as  the  FFT.  This  method  is  tested  on  three  images  demonstrating  various  compression 
ratios  currently  attainable. 


2  Derivation  of  Hyperdistribution  Transforms 


The  hyperdistribution  expansion  is  formulated  as  a  transform  which  is  similar  to  other 
orthogonal  transform  methods  with  the  addition  of  an  adjustable  parameter  which  controls 
the  shapes  of  the  basis  functions.  Begin  with  the  truncated  Rodriguez  hyperdistribution 
expansion  in  two  dimensions  as  derived  previously.  The  expansion  for  the  approximation  to 
a  function  f{x,y)  is  defined  as 
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with  coefficients 

( _ 1  \n+m  \Ti+Tn  *4-00  f+00 

L  L  «n(x/A)ff4WA)/(.,»)^x,i„.  (4) 

Insert  equation  (4)  into  (3);  cancel  (— separate  1/(2"^'”  n!m!)  into  square  roots, 
and  combine  1/A  with  each  differential.  Then  separate  the  equation  again  to  get  the  expan> 
sion 
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Now  substitute  for  the  function  f{x,y) 
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The  parameter  a  has  been  introduced  into  the  formulation  which  is  used  to  vary  the  shape 
of  the  wavelets  that  are  used  to  determine  the  approximate  function  (i.e.  reconstruct  the 
image).  The  shape  of  the  wavelets  will  affect  the  rate  and  type  of  convergence  of  the  expan¬ 
sion  series.  In  previous  Hyperdistribution  formulation,  a  can  be  considered  equal  to  one.  In 
the  case  were  a  is  equal  to  two,  the  wavelets  form  orthonormal  bases.  This  case  is  explored 
in  the  following  work 

Notice  that  1/A  appears  with  x  and  y  every  where  except  in  the  function.  The  variable 
substitution 


7 


X  =  Xu 
y  =  Xv 

dx  =  Xdu 
dy  —  Xdv 


produces  the  hyperdistribution  parametric  wavelet  transforms  with  the  expansion 

fe-i  fc-i  fj  /  \  IT  (  \  r --(u’+u*)! 
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on 
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3  Discrete  Formulation  of  Transforms 

The  equations  (12)  and  (13)  deal  with  the  expansion  for  a  continous  function,  but  we  are 
interested  in  discretely  sampled  functions.  An  actual  image  is  given  as  a  matrix  of  values 
which  can  be  considered  discrete  samples  of  a  function.  Consider  a  p  x  p  sampled  image 
where  i,j,  are  integers  and  Si,9j,  define  a  symmetric  coordinate  system  for  the  image,  define 

U  =  Sj  t  =  1 . . .  p 
V  =  Sj  j  —  \  ...p 

and  the  function  t  in  terms  of  these  integer  values 

tij  =  g{Xu,Xv)  g{x,y).  (15) 


From  equations  (12)  and  (13)  with  a  =  2  the  expansion  for  approximation  to  t  using  or¬ 
thonormal  wavelets  is 


fc-i  fc-i 
n"0  m=0 


(16) 


with  coefficients  r„„  =  using  the  discrete  version  of  equation  (13) 
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We  move  da  to  make  equations  (16)  and  (17)  symmetric  and  get  the  expansion 
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4  Matrix  Formulation  of  Transforms 

The  discretized  formulation  can  be  represented  in  more  conventional  form  by  using  a  matrix 
representation  similar  to  other  orthogonal  transforms.  Equations  (18)  and  (19)  suggest  using 
matrix  operations  for  doing  the  calculations. 

Consider  the  image  p  x  p  matrix 
<  tn  ti2  •  •  •  iip  \ 

<21  <22  •••  <2p  (20) 
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If  we  define  s.  k  x  p  Hermite  matrix 
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The  HD  coefficients,  r„„,,  can  now  be  caJculated  by  viewing  equation  (19)  as  the  following 
operation 


A  =  H  •  I  •  [Hp 


(22) 


which  creates  a.  k  x  k  IID  coefficient  matrix 


(Too  t-qi  ...  ro(fe_i) 

TlO  Tn  To(fc_ij 
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(23) 


We  can  reconstruct  the  image  by  viewing  equation  (18)  as  the  following  operation 


i  =  [H]^- A-H 


(24) 


which  creates  a  p  x  p  reconstructed  image  matrix 
^  tit  ti2  ■  •  •  tip  ^ 

A  *  A 

i  _  t2i  t22  "  •  t2p 

«  «  «  • 

«  •  «  • 

•  •  «  • 

A  A  A 
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(25) 


5  Implementation  Scheme 

An  algorithm  for  implementing  the  previous  formulas  to  images  is  developed  in  this  section. 
The  image  is  given  as  a  p  x  p  matrix  of  one  byte  integers. 

A  coordinate  system  is  determined  for  the  image  through  the  definitions 


3i  =  —L  +  (*  —  1)  da 

(26) 

(27) 

p-  1 

where  we  consider 

lij  =  I(s,-,Si) 

(28) 
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The  remaining  variables  to  be  determined  are  the  length  of  the  coordinate  space  of  the 
image,  (2  X  Z/),  and  the  dimensions  of  the  HD  coefficient  matrix,  [k  x  k).  This  gives  us  the 
parameters  L  and  k  to  adjust. 

We  predict  that  the  optimal  value  for  L  for  reconstruction  will  be  approximately  equal 
to  the  zero  closest  to  infinity  of  the  highest  order  wavelet  used  in  the  reconstruction.  This 
criterion  insures  that  all  the  wavelets  completely  fit  inside  the  coordinate  space  of  the  image 
and  that  wavelets  can  have  sign  variations  at  the  edge  of  the  coordinate  space  of  the  image. 
This  establishes  T  as  a  function  of  k.  Results  of  reconstruction  with  various  values  for  L 
substantiate  this  function  as  good  preliminary  criterion. 

The  zeros  of  the  nth  wavelet  are  equivalent  to  the  zeros  of  the  nth  order  Ilermite  poly¬ 
nomial,  since  this  is  the  only  part  of  the  equation  for  the  wavelet  that  can  change  sign.  For 
a  particular  choice  of  k  we  determine  L  by  finding  the  first  zero  from  infinity  of  the  {k  —  l)th 
order  Hermite  polynomial. 

The  compression  ratio  is  defined 


number  of  bits  in  the  image  matrix 

O ^  —  I  .i^  1^..,.— .  I  . . .  ,,  ^,1 . 

number  of  bits  in  HD  coefficient  matrix 
The  mean  squared  error  is  defined 
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MSB 


P  X  P  i^jTi 


p  X  p  (8bits)  p* 
k  X  k  (32bits)  4k^ 


(29) 


(30) 


8  Results  of  Compression/Decompressions 

The  first  image  analyzed  is  a  64  x  64  artificial  terrain  called  “Island”.  This  image  w.as  created 
by  combining  three  gaussian  functions  of  various  widths  and  strengths.  Gaussian  signals  are 
a  natural  test  case  for  reconstruction.  The  original  image  is  shown  by  surface  and  contour 
plots  in  figure  (1).  The  image  reconstructed  from  a  32  x  32  HD  coefficient  matrix  is  shown  in 
figure  (2).  The  image  reconstructed  from  a  16  x  16  HD  coefficient  matrix  is  shown  in  figure 
(3). 
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The  second  image  analyzed  is  a  256  x  256  image  called  “Gordon”.  The  original  image  is 
shown  in  figure  (4).  The  reconstructed  images  representing  compression  ratios  of  1,  4,  and 
16  are  shown  in  figures  (5),  (6),  and  (7). 

To  demonstrate  the  effects  of  the  parameter  L  we  analyze  the  reconstruction  of  “Gor¬ 
don”  with  higher  and  lower  then  predicted  values  for  L.  The  compression  ratio  for  each 
reconstruction  is  4.  The  predicted  value  for  L  was  used  in  figure(6).  Figure  (8)  shows  a 
reconstruction  with  L  one  half  the  predicted  value.  Figure  (9)  shows  a  reconstruction  with 
L  one  and  a  half  the  predicted  value. 

The  third  image  analyzed  is  a  512  x  512  image  called  “Liberty”.  The  original  image  is 
shown  in  figure  (10).  The  reconstructed  images  representing  compression  ratios  of  1,  4,  and 
16  are  shown  in  figures  (11)  and  (12). 

The  results  obtained  from  these  three  images  should  be  considered  in  the  nature  of  an 
existence  proof.  We  have  demonstrated  that  compressions  of  image  data  can  be  achieved 
with  HDT’s  and  recognizable  results  obtained  for  the  reconstructed  images.  It  is  particularly 
remarkable  that  these  results  have  been  obtained  for  unsegmented  images.  Additional  work 
to  be  undertaken  during  the  Phase  II  research  effort  will  improve  the  compression  ratio  and 
the  fidelity  of  reconstruction.  For  example,  conventional  techniques  such  as  image  segmenta¬ 
tion  and  image  processing,  combined  with  optimization  of  the  HD  wavelet  parameters,  can 
be  confidently  predicted  to  yield  continued  improvements  in  the  performance  of  the  II DT 
compression  algorithm. 


Figure  Captions 

•  Figure  1  Original  Image  “Island”  of  64  X  64  pixels;  (a)  surface  plot  (b)  contour  plot. 

•  Figure  2  Reconstructed  Image  “Island”  from  32  x  32  HD  coefficient  matrix  with  CR 
=  1  and  MSE  =  0.00451  .  (a)  surface  plot  (b)  contour  plot. 

•  Figure  8  Reconstructed  Image  “Island”  from  16  x  16  HD  coefficient  matrix  with  CR 
=  4  and  MSE  =  0.232  .  (a)  surface  plot  (b)  contour  plot. 
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•  Figure  4  Original  Image  “Gordon”  of  256  x  256  pixels. 

•  Figure  5  Reconstructed  Image  “Gordon”  with  CR  ==  1  and  MSE  “  13.6  . 

•  Figure  6  Reconstructed  Image  “Gordon”  with  CR  =  4  and  MSE  —  25.7  . 

•  Figure  7  Reconstructed  Image  “Gordon”  with  CR  =  16  and  MSE  =  59.7  . 

•  Figure  8  Reconstructed  Image  “Gordon”  using  low  L  with  CR  =  4  and  MSE  =  61.8  . 

•  Figure  9  Reconstructed  Image  “Gordon”  using  high  L  with  CR  =  4  and  MSE  —  1104 

•  Figure  10  Original  Image  “Liberty”  of  512  x  512  pixels. 

•  Figure  11  Reconstructed  Image  “Liberty”  with  CR  =  4  and  MSE  =  227  . 

•  Figure  12  Reconstructed  Image  “Liberty”  with  CR  =  16  and  MSE  =  476  . 
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Figure  3  Rcconslnu'lcd  Tnia^c  “Fsliuid”  from  10  x  10  III)  roi*iru  *u*nt  tiutlrix  with  (.U  -1  .mn  MSI 

=  0.232  .  (a)  surface  plot  (b)  contour  plot. 
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Figure  3b 
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Figure  7 
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7  Conclusions 


We  have  established  a  transform  description  of  the  hyperdistribution  approach  to  image 
processing,  the  hyperdistribution  transform  (HDT).  The  properties  of  the  HDT  have  been 
outlined  in  a  fashion  analogous  to  more  conventional  transforms,  such  as  the  Fast  Fourier 
Transform.  We  have  demonstrated  compression  of  both  synthetic  and  natural  images  using 
a  truncated  HDT  expansion.  For  unsegmented  natural  images,  compression  ratios  of  4:1  and 
16:1  were  demonstrated.  Standard  techniques  used  for  conventional  image  compression  al¬ 
gorithms  should  allow  further  improvements  in  the  performance  of  HDT  image  compression. 
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A  Computer  Code 


The  computer  program  used  to  perform  the  analysis  is  included  here.  The  language  used  is 
standard  ANSI  C.  The  utility  routine  used  are  from  Numerical  Recipes  Press  et.  al.  and 
are  not  included  in  this  report. 


/**♦♦*♦♦*♦♦**♦♦*♦♦*♦♦♦♦♦*♦***♦****♦♦**♦*♦♦♦**♦♦**♦****♦***♦♦♦***♦♦♦!*♦* 

*  This  program  is  used  to  calculate  the  Hyperdisdribution  transfrom  * 

*  of  an  image  and  to  reconctruct  the  approximate  to  the  image  from  * 


*  the  Hyperdistribution  transfrom.  • 

***  * 

♦  by  Thomas  W.  Drueding  date  :  3/9/91  * 

Aero-Mech  Dept .  * 

*  Boston  University  * 

♦  (617)  353-5260  ♦ 


#include "malloc . h” 

#include"nrutil . h" 

#include<stdio.h> 

#include<niath .  h> 

#defino  PI  3.14159265359 

void  image_read() ,image_write() ,init_grid() ,make_hdimg() ,make_rlimg() ; 
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main (argc , nrgv ) 
int  argc; 
char  ♦argv[]  ; 

■C 

float  ♦♦imaga, **recon, **hdiing,**trans, limit ; 
int  p.lc; 

/*  Inititialize  paramtera  ♦/ 
if  (argc  <  3)  { 

printf ("Usage:  hd  (p)  (k)  (limit)  \n"); 
oxitO  j 

} 

Bscanf  (argvCl]  ,"y.d",tp) ;  /♦  x/y  size  of  real  image  (assume  square  image)  */ 

sscanf  (argvC2]  ,"Xd",ft]c)  :  /♦  x/y  size  of  HD  image  (assxime  square  image)  */ 

sacanf (eo-gvCS] ."Xf " .klimit) ;  /♦  edge  of  coordinate  system  on  the  image  */ 


/*  Allocate  dynamic  memory  for  matrices  */ 

image  =  matrix(l,p,l,p) ;  /♦  image  :  original  image  */ 
recon  =  matrix(l,p,l ,p) ;  /♦  recon  :  reconstructed  image  */ 
hdimg  =  matrix(0,k-l,0,k-l) ;  /♦  hdimg  :  HD  image  (HD  transform  of  imago)  */ 
trans  =  matrix(0,k-l,l,p) ;  /*  trans  :  HD  transformation  matrix  */ 

/*  PROCESS  */ 

image.read (image, p," image. b" ) ;  /*  road  original  image  from  file  */ 
init_grid(trans, k,p, limit) ;  /*  setup  HD  Transformation  matrix  ♦/ 
make_hdimg (image, hdimg, trans, k, p) ;  /*  calculate  HD  transform  of  image  */ 


mako_rlimg(rocon, hdimg, trans, k, p) ;  /♦  calculate  image  from  HD  transfrom  */ 
imago_writo(rocon,p, "rocon.b") ;  /♦  write  reconstructed  image  to  file  */ 

> 
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/*  Read  BINARY  image  from  file  */ 
void 

image_read( image,  size,  filename) 
float  **image; 
int  size; 
char  ♦filename; 

■c 

FILE  ♦fp; 
int  i,j; 

if  ((fp  =  f  open  (filename,  "rb”))  ==  NULL)-C 
printf ("Cannot  open  file  y,a\n",  filename); 
exit(O) ; 

> 

printf ("NnReading  y.8...\n",  filename); 
for  (i=l;  i<=8ize;  ++i) 
for  (j=l;  j<=8ize;  ++j) 

image [ij[j]  =  ((float)  getc(fp))  -  127.0 
f close (fp) ; 
printf (". . .done\n") ; 
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/♦  Write  BINARY  image  to  file  */ 
void 

image_write( image,  size,  filename) 
float  ** image; 
int  size; 
char  *filename; 

i 

FILE  *fp; 
int  i,j; 
float  val; 

if  ((fp  =  f  open  (filename,  "wb"))  ==  HULL)-C 
printf( "Cannot  open  file  XsW,  filename); 
exit(O) ; 

} 

printf("\nWriting  y.s..An",  filename); 
for  (i=l;  i<=8ize;  ++i) 
for  (j=l;  j<=size;  ++j){ 
val  =  image  [i3tj]  +  127.0; 

^f  (val  <  0.0)  val  =  0.0;  /*  check  bounds  ♦/ 

else  if  (val  >  255.0)  val  =  255.0; 

putc(  (char)  val  ,  fp) ;  /•  convert  to  binary  and  write  */ 

> 

fclo8e(fp) ; 

pr intf ( " . , , done\n" ) ; 
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/*  setup  HD  Transformation  matrix  */ 
void  init_grid(trans ,k,p, limit) 
float  ♦♦trans, limit ; 
int  k,p; 

■C 

int  n,i; 

float  8,ds,Const,Tn2,Tnl,TnO; 

printf ("\nlnitializing  grid  ..An"); 
ds  =  2.0  ♦  limit  /  (p-1.0); 

Const  =  aqrt(  ds  /  sqrt(PI)  ); 
for  (i=l;  i<=p;  i++){ 

8  ®  -limit  +  (i-1)  *  ds; 

Tnl  *  0.0; 

transCOJCi]  *  TnO  =  expC-s  *  8/2)  ♦  Const; 
for  (n=l;  n<=k-l;  n++  ){ 

Tn2  *  Tnl; 

Tnl  »  TnO; 

trana[n][i3  =  TnO  =  s  ♦  sqrtC (float )2/n)*  Tnl  -  sqrt ( (float) (n-l)/n)  *  Tn2 

} 

> 

printf ( " . . . done\n" ) ; 

} 
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/*  calculat«  HD  transform  of  image  ♦/ 
void  make.hdimgCrlimg , hdimg , trans , k, p) 
float  ♦*rlinig,**hdimg,**tran8; 
int  k,p; 

i 

int  i,j,n,m; 

printf ("VnMaking  HD  image. . An") ; 
for  (n=0;  n<=k-l;  n++) 
for  (m=0;  m<=k-l;  m++) 
for  (i=l;  i<=p;  i++) 
for  (j=l;  j<=p;  j++) 

hdis^Cn]  [m]+=rliing[i]  [j]*tran8[n]  [iJ^tranaCm]  [j]  ; 
pr intf ( " . . . done\n" ) ; 

> 

/♦  calculate  image  from  HD  transfrom  ♦/ 
void  make _rlimg(rlimg, hdimg, trans, k,p) 
float  ♦*rlimg, ♦♦hdimg, ♦♦trans; 
int  k,p; 

{ 

int  i,j,n,m; 

printf ("\nMaking  RL  imago ... \n" ) ; 
for  (i=l;  i<=p;  i++) 
for  (j=l;  j<=p;  j++) 
for  (n=0;  n<=k-l;  n++) 
for  (m=0;  m<=k-l;  m++) 

rlimgCi] tj]+=hdimg[n] [m] ♦trans [n] [i] ♦trans [m] [j]  ; 
printf (". . .dono\n"); 

} 


35 


B  Preprints  of  Publications 

la  this  appendix  we  attach  two  preprints  of  publications  submitted  under  this  contract. 


Hyperdistributions  I:  One  Dimensional  Analysis 

T.  Drueding^,  T.  Fang^,  R.  Hohlfeld*,  W.  Ma*,  G.  Sandri* 

*  Boston  University  College  of  Engineering 
no  Cummingtion  Street,  Boston,  MA  02215 

*  Massachusetts  Technological  Laboratory,  Inc. 

330  Pleasant  Street,  Belmont,  MA  02178 


Submitted  for  Publication:  SIAM  Journal  for  Numerical  Analysis 


Hyperdistributions  I:  One  Dimensional 

Analysis 


We  develop  in  tiiis  paper  and  the  following  paper  a  new  technique  for  the 
calculation  of  convolution  products  and  their  inverses.  This  is  accomplished 
by  constructitig  a  class  of  singular  “functions”,  hypcrdistrihution,  that  form 
a  closed  algeiiraic  field  with  the  convolution  product  as  the  multiplicative 
operation.  In  this  paper  we  consider  functions  of  a  single  variable.  This 
one-dimensional  algebra  can  be  applied  to  signal  processing.  Furthermore, 
in  this  paper  we  use  the  construction  of  hyperdistributions  to  obtain  a  novel 
parametric  approximation  method.  We  demonstrate  the  use  of  our  approxi¬ 
mation  method  with  simple  examples. 
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Introduction 


In  the  sequel,  we  develop  a  novel  tool  ff>r  the  ralnilatifui  of  eoiivolution 
produrls  and  their  inverses;  hyperdistrihiitions  which  give  simple  applica 
tions  to  the  two  related  areas  of  signal  analysis  and  image  processing. 

Hyperdistrihiitions  are  singular  “functions’*,  defined  and  constructed  he 
low,  that  we  find  to  form  a  algebraic  ahelian  field  with  the  convolution  prod 
uct  as  the  multiplicative  operation. 

Signal  analysis  is  a  natural  application  of  hyperdistrihution  in  one  (one 
dimensional)  variable.  Image  processing  is  a  corresponding  application  for 
hyperdistrihiitions  in  a  two  dimensional  variable.  Tomography  corresponds 
to  three  dimensional  variable  and  the  budding  field  of  space-time  processing 
corresponds  to  four  dimensions. 

The  outline  of  the  paper  is  as  follows.  In  section  I  we  introduce  hy¬ 
perdistributions  heuristically.  ive  discuss  the  convolution  group  and  derive 
a  remainder  theorem.  In  section  11  we  construct  rigorously  hyperdistribu- 
tions  by  a  modified  Hermite  polynomial  expansion  and  we  use  the  tools  of 
the  Christoffel-Darboux  theory  to  obtain  sufficient  conditions  for  conver¬ 
gence.  In  section  III  we  show  that  Clauss’  multipole  expansion  is  obtained 
explicitly  as  a  simple  application  of  the  hyperdistribiition  inverse.  Finally,  in 
section  IV  we  expand  a  gaussian  function  in  terms  of  derivatives  of  a  different 
gaiissiau  to  demonstrate  the  use  of  our  new  parametric  expansion  and  the 
concurrent  minimization  of  error.  The  logical  interconnection  of  the  sections 
is  shown  in  Fig  1. 

We  note  that,  in  effect  we  introduce  a  method  for  establishing  and  ap¬ 
proximating  solutions  of  integral  equat'ons  of  convolution  form  based  on  the 
use  of  this  new  class  of  singular  functions.  We  demonstrate  with  examples 
that  there  are  rases  for  which  our  method  is  applicable,  but  Fourier  trans¬ 
form  methods  fail.  Our  method  requires  the  cab  illation  of  the  moments  of 
the  given  functions  and  of  the  kernel  rather  than  calculation  of  their  Fourier 
coefficients.  This  property  motivates  consideration  of  examples  for  which  our 
method  is  clearly  preferable  to  Fourier  transform  techniques.  Applications 
are  given  with  emphasis  on  image  deconvolution  and  the  analysis  of  turbulent 
flows. 
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1  Heuristic  definition  of  Hyperdistribution. 

We  introduce  a  genera!  approximation  whose  integral  properties  are  the  focus 
of  interest.  Our  approximation  displays  in  configuration  spare  the  properties 
of  the  classical  moment  generating  expansion  for  the  Fourier  transform  of 
the  prohahility  distribution.  Moments,  and  even  shapes,  are  shown  to  be 
captured  well  by  our  expansion.  In  addition,  our  expansion  allows  consider 
ation  of  “functions”  which  are  more  singular  than  temperate  (  i.e.  Fourier 
analyzable)  distributions,  but  that  ran  be  represented  by  infinite  sequences  of 
distributions.  We  call  these  “functions”  generalised  distributions.  We  show 
that  generalised  distributions  form  the  appropriate  framework  for  carrying 
out  the  process  of  deconvolution,  in  fact,  we  give  a  st'^aightforward  algorithm, 
the  “Bochner  Martin  algebra”,  to  compute  explicitly  the  convolution  inverse 
of  any  generalised  distribution.  Applications  of  the  method  to  simple  opti¬ 
cal  deconvolutions  are  shown  to  be  straightforward.  These  applications  have 
been  motivated  originally  by  our  studies  of  optics  and  turbulent  flows.  W? 
have  since  become  aware  of  the  much  wider  applicability  of  our  methods. 

1.1  Taylor  and  moment  expansions 

The  theory  of  generalised  distributions  is  built  on  ideas  related  to  Dirac  's 
delta  function,  which  is  technically  a  “generalised  function”.  The  Dirac  delta 
function.  S(t).  has  the  properties 

I  fi{r)dx  -  ]  ( 1 ) 

*/  -  OC 

and 

/(.r)  -  f  S{t  x')f{T')dT'  (2) 

J  OO 

for  any  suitably  smooth  function,  /(tV  We  term  equation(2),  “Dirac's  iden 
tity”. 

The  Dirac  delta  function  is  symmetric  in  its  argument,  i  e.  f>[x  x') 
x)  and,  since  the  Dirac  delta  function  is  a  generalised  function  .  we 
may  laylor  expand  /'(r'  r)  about  x’. 

f>{x'  x)  -  x^’{x')  \  xH"{x')  I  ...  I  /?’’(t)  (3) 
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Here  6'(t)  “  d8(T)ldr,  8  "(t)  d‘8[T)/dT' .  elr.  l  iiis  expression  aiui  similar 

ones  below  are  assumed  lo  hold  under  intepralion.  F'qnation  (3)  allows  one 
1o  compute  a  “local  '  approximation  t<»  /(r),  since  if  we  stjbstitute  this  ex 
pansion  into  Dirac  s  identity,  we  recorer  the  usual  laylor  expansion  of  /(x) 
about  T  “  0. 

f(T)  -  f{Q)  t  A  T^/2!/"{n)  4  ...  {A) 

This  approximation  is  local  in  the  sense  that  it  requires  derivatives  of  /{t) 
at  a  single  point,  x  —  0,  and  in  general  has  a  limited  radius  of  convergence. 
On  the  other  hand,  if  we  expand  6(x  x')  about  x.  we  have 

^(x  -  x')  -  /)(x)  -  t'8'{t)  t  x'V2!^"{x)  ^  {•'’) 

When  this  series  is  substituted  into  Dirac  's  identity,  we  obtain 

fix)  -  4  M^/2U"{x)  4  ...  4  /?"(x)  (6) 

The  coefficients  A/",  defined  by 

yi/"  -  f  x^f{x)dx  (7) 

are  simply  the  moments  of  the  function  f{x).  Therefore  equation  (6)  is  an 
approximation  of  /(x)  involving  global  information  about  f{x),  that  is,  the 
moments  of  the  function. 

Th  is  then  may  be  taken  as  a  motivation  for  our  definition  of  a  generalised 
distribution  as  a  function  which  may  be  written  in  the  form, 

fix)  -  ^a„rix)  (8) 

n  ~  0 

The  values  are  coeffirients  given  by 

a„-i  iriir/n!  (9) 

Note  that  equation  (8)  is  equivalent  to  the  familiar  moment  gefierating  ex 
pansion  of  probability  theory  (see  equation  (6)). 
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V/  , 


1.2  The  Convolution  Group 

If  we  have  two  generalised  distributions  /,  and  /j,  their  linear  combination 
+  /'/i  ”  generalised  disiribniion  (where  A  and  //  are  real  coef¬ 

ficients).  The  pth  derivative  of  a  generalised  distribution,  (f'’/(x)/rfr’'  - 
V'’  /(x).  is  a  generalised  distribution.  Also,  the  convolution  of  two  gener¬ 
alised  distribution  /i  *  /2  is  a  generalised  distribution.  These  may  all  be 
thought  of  as  “clostire  properties"  of  generali.sed  distributions. 

Generalised  distributions  allow  us  to  make  an  effective  computation  of 
the  convolution  inverse.  Given  a  generalised  distribution  /,  the  desired  con 
volution  inverse  /n|/|  satifies 

(10) 

Here  8  represents  Dirac  s  delta  function,  which  is  the  identity  of  the  convolu 
tion  operation.  We  shall  show  by  our  construction  that  /n|/!  is  a  generalised 
distribution.  Writing  the  convolution  explicitly, 

/  *  Ini/!  -  r  dx'J{x')]n\f]{T  ■  x')  ^(x)  (11) 

J  -  rx) 

It  is  now  necessary  to  compute  the  product  of  the  sums  and  match  coeffi¬ 
cients.  Taking,  /(x)  -  o,,  V"  (^)’ ^nd  /n|/{x)]  -  V"  ^(^)i  see 

that  the  computation  of  the  convolution  inverse  is  effectively  the  determina¬ 
tion  of  a  collection  of  values,  given  a  set  of  o„  values.  Substituting  into 
equation  (11), 

/  rfx'(^o,v'"^(^'))(Z^V'’^(^  -=r'))-Ii{x)  (12) 

•'  p--0  <7-0 

Noting  that  r'  ~  Vr  ^nd  that  v'’^(^)*  V’^(^  ^')  “  V'’^’^(3‘)< 

(12)  may  finally  be  redticed  to 

P  7 

or  equivalently  with  r  p  i  q 

oc  no 

pv'^(^)-  Ht)  (11) 

p  Or  0 


Matching  coefficients  on  the  left  hand  and  right  hand  sides  implies  that  only 
the  T  0  term  survives.  The  result  is  a  linear  system  of  equations  for  the  h 
values  in  terms  of  the  n  values.  It  is  easiest  to  see  the  behavior  hy  writing 
the  first  few  equations  in  this  linear  system. 


flo^O  1 

flo^i  -1  fli^o  ~  0 

Q0^2  t  t  ^2^0  0 

t  0  ihj  -|  t  n^hf)  0 


(l.'i) 


and  so  forth.  Thus,  we  can  see  that,  feo  ~  1 /fln.  - 

02/flo.  etc.  This  completes  the  computation  of  the  terms  of  the  convolution 
inverse. 


1.3  Fourier  Transforms  and  simple  examples 

Lastly,  we  shall  want  to  consider  the  Fourier  transform  of  a  generalised  distri¬ 
bution.  Again  taking  /(r)  --  oO„  V"  ^(r),  we  can  immediately  evaluate 
the  Fourier  transform  as: 

J  infty 


where  for  the  last  step  we  have  used  —  i’^k.  The  Fourier  transform  of 

a  generalised  distribution  may  be  seen  to  be  a  power  series. 

The  Fourier  transform  of  a  function  is  the  “moment  generating  function" 
because  of  equation  (9).  Thus  the  basic  requirement  for  the  ^-alidity  of  writing 
a  function  as  a  generalised  distribution  is  that  its  moments  be  finite,  or  more 
stringently,  that  its  Fourier  transform  be  real  analytic. 

Generalisation  of  this  discussion  to  a  two-dimensional  generalised  distri 
bution.  as  would  be  required  for  most  types  of  imaging  data,  is  relatively 
straightforward.  This  generalisation  retains  the  mathematical  properties, 


^  a.,  /  V" 

n-O  * 

n^O  '  ^ 


(16) 


n- 0 
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particularly  for  convolutions  and  Fourier  transforms,  that  make  generalised 
distributions  attractive  for  signal  processing. 

Simple  explicit  examples  of  convolution  inverses  are  obtained  from  the 
standard  Creen  "s  functions.  The  one-dimensional  Helmholtz  equation,  for 
example,  can  be  treated  equally  easily  with  generalised  distributions  and  with 
Fourier  transforms  to  obtain  both  the  Creen ’s  functions  and  its  convolution 
inverse.  Thus,  we  have 


{d'^ldx'^  -4  ]/4)e''''/^  rr  -8[t)  (17) 

Application  of  equation  (15)  shows  that  only  bo  and  hj  contribute  to  the 
convolution  inverse  and,  in  fact,  W'riting  /  —  lnv{f], 

=  l/26(r)  -  6"(r)  (1^^) 

This  result  is  easily  reproduced  by  using  Fourier  transforms  and  is  easily 
verified  using  equation  (10). 

In  contrast,  the  convolution  inverse  of  the  Gaussian  function  cannot  be 
obtained  ar  a  Fourier  transform,  but  it  is  easily  obtained  using  generalised 
distributions.  It  may  be  verified  that  for  a  Gaussian  function  of  width  <r,  the 
convolution  inverse  is  readily  obtained  from  equation  (15),  while  calculation 
using  Fourier  transforms  leads  to  a  divergent  result.  Calculation  of  moments 
and  use  of  equation  (15)  yields 

a2n  =  {<r/\/2)"/n! 

b2n  -  (  l)"(o/v^)"/n!  (19) 
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Rigorous  Construction  of  Hyperdistribu¬ 
tions 


In  this  section  we  first  modify  the  classical  expansion  in  the  Ilermite  polymi- 
als.  This  expansion  is  “superior”  to  a  power  series  expansion  in  that  the  terms 
are  orthogonal,  making  the  error  orthogonal  to  the  approximation.  Our  mod¬ 
ification  maintains  this  advantage.  Fnrlhurmore.  the  standard  Christoffel- 
Darboux  analysis  gives  a  very  useful  sufficient  condition  for  convergence  that 
our  expansion  inherits  from  the  Hermite  polymial  expansion.  Since  our  ex¬ 
pansion  utilises  the  Rodriguez  formula  for  the  Hermite  polymials,  we  call  it 
the  “Rodriguez  expansion”. 

2.1  Modified  Hermite  Expansion:  Rodriguez  expan¬ 
sion 

In  effect,  we  discuss  here  a  systematic  pointwise  approximation  of  generalised 
distributions.  This  approximation  method  is  analogous  to  the  method  de¬ 
veloped  by  Temple  to  approximate  distributions  by  smooth  functions.  The 
Temple  method  has  become  widely  known  through  Lighthill  's  monograph: 
Fourier  Series  and  Generalised  Functions. 

Generalised  distributions  are  approximate  by  a  modification  of  the  clas¬ 
sical  expansion  in  Hermite  polymials  (hereafter  called  the  “Hermite  expan¬ 
sion”)  w'hich  is  suggested  by  the  Rodriguez  formula: 

(20) 

where  V  denotes  the  derivative.  For  simplicity,  we  consider  here  the  one¬ 
dimensional  problem,  so  that  sj  is  the  ordinary  derivative  with  respect  to  the 
variable  x.  To  see  how  the  Hermite  expansion  can  be  transformed  into  an 
approximation  for  generalised  distributions,  consider  a  function  /(r)  which 
we  want  to  represent  as  a  generalised  distribution.  Multipy  /{x)  by 

then  expand  the  resulting  expression  in  terms  of  the  "scaled  Hermite  polymi¬ 
als 

//„^(t)  a  -rL{x/\)  (22) 
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whose  definition  is  justified  by  the  formula  (25)  below.  VVe  obtain  the  ex 
pression 

n-0 

Multiply  both  sides  of  the  equation  by  the  normalised  CJaussian  introdured 
by  Temple,  i.e. 

(24) 

And,  observe  that  by  rescaling  the  Rodriguez  formula  (20)  we  ran  write 

(-l)"/f„V)^>(a-):.  (2fi) 

The  final  result  is 

/(»)  "  £  Ci  V"  <>(*)  (26) 

n— 0 

where  we  have  introduced  the  coefficients 

“2  =  (22) 

^  Tl  •%/■*“  oo 

When  we  let  A  tend  to  zero,  we  have  a  representation  of  f{x)  as  a  generalised 
distribution.  We  call  the  expansion  (26)  the  ’’Rodriguez  expansion  for  f{x) 
wuth  width  A”.  The  width  parameter  is  a  novel  feature  of  the  Rodriguez 
expansion  when  compared  to  standard  expansions  in  complete  sets  of  basis 
functions.  The  standard  expansions  do  not  contain  a  free  parameter.  The 
advantages  of  the  Rodriguez  expansion  will  be  demonstrated  below  in  the 
context  of  the  theory  of  generalised  distributions  and  of  convolution  inverses. 

We  consider  as  an  example  a  function  familiar  from  the  analysis  of  tur¬ 
bulence  spectra,  the  “Ogura”  function 

(28) 


With  N  determined  by  the  normalisation  condition, 


we  have 


-  (l/2)r(l/4)  1.812R  (30) 

A  number  of  graphs  are  now  used  to  illus1ra<e  (he  three  main  advantages  of 
the  Rodrigue/  expansion. 

1.  The  Rodriguez  expansion  converges  point  wise.  (See  Fig  3.)  The  fler- 
mite  expansion  also  converges  pointwise,  but  the  convergence  is  slow 
for  large  values  of  the  argument.  In  addition,  the  Hermite  expansion 
exhibits  "  whipping  tails”  for  large  values  of  the  argument. (See  Fig  6.) 

2.  The  pointwise  convergence  of  the  Rodriguez  expansion  is  independent 
of  A,  but  the  rate  of  convergence  can  be  optimised  by  a  proper  choice 
of  A.  (See  Fig  3) 

3.  In  addition  to  pointwise  (i.e.  local)  convergence,  the  Rodriguez  expan¬ 
sion  exihibits  a  global  convergence  manifested  by  accurate  represen¬ 
tation  of  moments  of  the  function  being  approximated.  To  illustrate 
this  global  convergence,  since  the  function  we  are  considering  here  is 
symmetric  in  we  have  found  it  convenient  to  consider  ’’one-sided 
moment”  plots  to  exihibit  the  convergence  of  our  generalised  distri¬ 
bution  approximations.  These  are  integrals  over  positive  x  only,  i.e. 
p„(x)  =  /,J'(x')”/(x')<ix'.  The  one-sided  moments  of  order  n  are  given 
with  excellent  accuracy  by  numerical  integration  of  the  approximation 
with  n  -f  1  terms.  (See  Figs  i  and  5.)  Again  in  sharj-  contrast,  the 
Hermite  expansion  gives  divergent  moments.  (See  Fig  7)  Note  that  the 
horizontal  axis  in  the  one-sided  moment  plots  is  the  argument  x. 

To  summarise,  the  Hermite  expansion  exhibits  pointwise  convergence,  albeit 
with  “whipping  tails”.  By  contrast,  the  Rodriguez  expansion  exhibits  point- 
wise  convergence  without  “whipping  tails”.  As  a  consequence,  the  Rodriguez 
expansion  also  yields  a  global  approximation  by  accurately  representing  the 
moments.  Furthurmore,  the  Rodriguez  approximations  has  an  adjustable 
parameter  w'hich  may  be  selected  to  optimise  the  rate  of  convergence.  The 
Christoffel- Darboux  theory  gives  the  necessary  condition  for  convergence  in 
appropriate  L^.  For  the  Hermite  expansion  we  have 


•V  • 


Inserting  the  transformation  lhat  leads  to  the  Rodriguez  expansion,  we  find 


J  -  nc) 


dr  < 


oc 


for  the  corresponding  definition  for  the  Rodriguez  series. 


(32) 


2.2  Definition  of  Hyperdistributions 

In  order  to  describe  the  process  of  antidiffusion,  we  have  introduced  a  class 
of  highly  singular  “functions”,  that  are  precisely  defined  in  this  section.  Our 
process  for  defining  hyperdistributions  parallels  the  Temple  definition  (  gen¬ 
eralised  function)  as  a  good  sequence  of  good  functions.  Good  functions  are 
smooth  and  tapered.  More  precisely,  they  are  total  point  functions  that  are 
differentiable  to  all  orders  (G**),  and  decay  at  Too  faster  than  any  powrer. 
Good  functions  play  the  role  of  “testing”  a  sequence  of  good  functions  for 
weak  convergence.  In  fact,  a  sequence  of  good  functions  is  a  distribution  if 

lim  /  <^(r)/n(x)c/x  <  oo  (33) 

J  -oc 

for  all  good  functions  d». 

Since  we  conceive  of  hyperdistributions  as  “generalised”  distributions,  we 
are  in  fact  implementing  a  second  order  generalisation  of  functions.  Con¬ 
sequently.  we  need  a  double  test  as  a  convergence  criterion.  We  implement 
this  criterion  by  introducing  very  good  functions  Ga(x)  with  the  following 
properties: 

1.  Ci^ix)  is  smooth,  that  is,  differentiable  to  all  orders,  C^. 

2.  Ga(x)  is  essentially  compact,  i.e.  it  has  a  gaussian  decay  at  -fee  : 
Ga(x)  - 

W'e  will  assume  for  convenience  that  Ga  is  normalised  to  unity: 

Gn{T)dx  -  1 

W'e  define  the  width  of  G\  by 

/♦  OC 

(x  t)"G  f^{T)dT 

no 


(3-1) 

(3r,) 


II 


A  primary  example  of  very  good  functions  is  the  gaussian,  which  we  denote 
by  6x{x): 

We  now  introduce  a  sequence  of  very  good  functions  defined  by 

h{x)  (37) 

fc-O 

The  sequence  where  A  is  a  nonnegative  real  and  n  is  a  natural 

number,  is  a  good  sequence  if,  for  all  good  functions  (f>  and  for  all  very  good 

functions  G\.  there  exists  a  Aq  such  that,  for  all  A  >  Ao, 

lim  r*  (f>{x){n^  *  C?A  )(x)d3-  <  oc  (38) 

r.0-^ - 

We  note  that  6{x)  are  hyperdistributions.  The  sum  (hyperdistribution) 

J;aAV''^(3-)  (39) 

Jb-O 

can  thus  be  viewed  either  as  a  sequence  of  “good”  distributions  as  A  -»  0: 

(40) 

k-o 

or,  as  n  “  »  oo,  as  sequence  of  good  functions: 

^  V*"  ('ll) 

k-0 

The  latter  representation  is  a  Rodriguez  expansion.  The  Rodriquez  formula 
for  the  Hermite  polymials  can  be  used  to  show'  the  derivatives  of  a  gaussian 
form  a  complete  set  of  orthogonal  polynomials  in  an  space.  And  thus 
the  Rodriguez  expansion  yields  a  very  useful  point  function  approximation 
to  any  hyperdistribution: 

(12) 

*-0 

where  ffnix)  denotes  the  Ilermite  Polynomial  in  x  of  order  n. 
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3  The  Multipole  Expansion 

We  start  with  the  familiar  Poisson  equation  of  potential  theory. 

(43) 

This  equation  is  rewritten  with  the  help  of  the  (infinite  domain)  Creen  ’s 
function 

C{t)  rr  _-l,  ^  ^  W 

47rr 

VVe  can  then  rewrte  the  “potential’’,  4>.  in  terms  of  the  “charge  distrihotion”, 


Introduce  Q  with  the  property 


(f>  -  G  *  p 


Q  *  p  -  S 


Convolving  hoth  sides  of  eq(45)  with  Q  and  using  the  commutative  and 
associative  properties  of  the  *  product,  we  find 

Q  *  <f)  —  Q  *  (G  *  p) 

=  (Q*p)*(^ 

=  G  (47) 

Solve  eq(13)  for  0  in  terms  of  the  given  G  by  computing  the  convolution 
inverse, 

4>{t)^.^\^S7’’G{x)  (48) 

fc-o 

which  is  Gauss’  multipole  series  with  coefficients 


At  =  j  x^''p{T)d\ 


where  (gi  denotes  tensor  product  and  a-**'  is  the  tensor  powder  of  the  3- vector 


Substituting  (49)  into  (48)  we  find  a  familiar  expression 

f  1 1*'  /•  - 1 

<^(®)  ^  ~  Y.  ^  li  ( /  piy)<^^y)  ® 

which  is  a  standard  result  in  potential  theory. 


(50) 


We  can  interpret  the  Rodriguez  expansion  as  a  generalised  mnitipole  ex¬ 
pansion  that  includes  a  size  or  radius  parameter,  A.  Thus  our  mnnopole 
generates  the  point  source,  which  is  a  Dirac  delta  function  as  a  gaussian  of 
width  A.  We  recover  the  standard  multipoles  when  A  ►  0.  Our  expansion 
has  the  form 


f(x)  rr.  a^f^xir)  )  V  /*a(3')  ‘  flj  ^a(3-)  4  ... 

(51) 

The  ingredients  of  the  expansion  are  the  basis  functions 

(52) 

and  the  coefficients 

Z,-'  z.^ 

Of,,  o, ,  02, ... 

(53) 

all  of 

whicli  depend  on  the  size  parameter  A.  The  basis  functions 

are,  explic- 

itly 

Sx{x)  =.  - 

\/irX 

(54) 

VfA(T)  -  -2x~—r- 
y/nX 

(55) 

V“  «,(»■)  -  (iz-V*"  -  2/.'’)^^ 

v/ttA 

(56) 

These  functions  are  plotted  in  Figs  8,9  and  10  for  3  values  of  A  which  show 
how  the  basis  functions  tend,  as  A  0,  to  a  monopole,  a  dipole  and  a 
quadrupole  respectively.  The  corresponding  expansion  coefficients  are  given 

Fy 

,  „  ^  4  oo 

4  -  /  f{^)dT  (57) 

J  —  no 

This  coefficient  represents  the  total  area,  or  alternatively,  the  total  mass  or 
charge  of  the  source  function.  Also 

-  f  xf[x)dx  (58) 

J  OO 

which  represents  the  net  dipole  of  the  source.  In  mechanical  terms  defines 
the  center  of  mass  of  the  source  function.  Finally 

a^^ff'x^f{x)dx  (59) 

•f  -IX* 


M 


which  represents  the  quadnipole  of  the  sotirce  distributions, 
terms  this  is  the  moment  of  inertia  of  the  mass  distribution, 
basis  function  and  coefficients  have  analogous  interpretations. 


In  mechanical 
Higher  order 
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4  Expansion  of  a  Function  in  Rodriguez  Se¬ 
ries  and  its  Optimization. 

Oiir  expansion  preserves  the  riassiral  prf)periirs  ifiat  are  (ierivablc  for  orthog¬ 
onal  polynomial  exapnsions  (as  opposed  to  power  series),  bnl  also  adds  an 
important  new  feature;  the  size  or  radius  parameter  that  generates  the  Gaus 
sian  picture  of  point  monopole,  dipole,  quadrupole.  etc,  to  a  scenario  in  which 
we  can  allow  for  extended  sources.  We  find  from  the  convergence  condition 
derived  from  the  Christoffel  Darboux  thoery  in  section  II  that  convergence 
holds  generally  on  a  semi-infinite  range  of  A.  This  remarkable  freedom  is  ex¬ 
ploited  in  this  section  to  optimise  the  rate  of  convergence  of  the  expansion. 
This  optimisation  results  in  determining  the  value  of  A  for  which  a  minimum 
number  of  terms  is  determin-  d  in  order  to  obtain  a  given  tolerance.  We  mea¬ 
sure  the  tolerance,  as  usual,  by  the  least  square  fit  integrated  over  the  entire 
function.  We  define  the  error 

t{N,X)  ^  r'drU^T)  -  i;  oi  V"  Ml))’  (60) 

■'-«  n-O 

We  apply  the  formula  to  a  standard  gaussian,  i.e.  we  take 

/(i)  =r  c  (61) 


we  then  find,  after  some  calculations 


.Vr-O 


(  l)*^2''  -=-"(2n  -  2k  -  1)!! 


-I  (  !)"2  ^"/n' 


k\{2n  —  2ky. 


where  the  double  factorial  contains  only  odd  terms.  We  can  now  express  the 
standard  gaussian  in  the  Rodriguez  form: 


.-'Vvv  -  E 


(>  -  ^’)" 


22"n! 


V  hir) 


A  proof  based  on  useful  linear  operator  relations  is  as  follows.  We  can  check 
by  differentiation  in  f  and  t  that 
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We  now  let 

4/  A^ 

(6.'-,) 

As  a  consequence,  w-e  can  write 

-  h^{T) 

(66) 

V^rA 

We  now  consider  the  identity 

r- 

(67) 

expanding  the  R.H.S. 

__  n  A^t"  _ 

~  c-'’/v^  -  Y.  v’"  «>(*)  (68) 

n-0  ^  ^• 


The  methf-H  of  proof  will  he  used  later.  Furlhurmore,  it  provides  a  welcome 
check  on  the  rather  difficult  calculations  of  coefficients.  We  see  hy  inspection 
that  A  —  1  ib  optimal.  Fig  1'  shows  the  trend  of  the  error  function  obtained 
numerically.  Fig  12  shows  the  error  function  on  an  expanded  scale.  The 
behavior  is  smooth  and  gradient  searches  promise  to  he  straightforward.  We 
also  consider  the  function 

fix)  -  co.o{Kx)^  (69) 

•v/ttA 

The  overall  behaviour  is  shown  in  Fig  13  and  an  enlargement  in  Fig  14.  The 
function  exhibits  rapid  variation.  Fig  7  shows  the  choice  of  parameters  K  -■  4 
and  —  3.  Fig  16  shows  the  behaviour  of  the  error  and  the  corresponding 
optimisation.  In  Fig  17  and  18  we  show?  the  result  of  deconvoluting  the  letter 
T  with  hyperdistribution  algebra  after  it  was  “smeared”  with  a  gaussian 
filter.  In  Fig  19,  we  simulate  the  crossberg  Todororich  neural  network  for 
the  cornsweet  effect  with  hyperdistributions. 
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5  Conclusions 


We  have  seen  that  generalisefl  disirihniion  prf)vide  a  method  for  solving 
Fredholm  integral  equations  of  the  ronvoliifion  type  which  is  competitive 
with  the  current  methods  that  employ  Fourier  transforms.  We  have  also  seen 
that  generalised  distributions  can  be  approximated  numerically  by  sequences 
of  smooth  functions.  This  procedure  is  analogous  to  that  of  approximating 
Dirac  f>  functions  by  sequences  of  narrowing  Gaussian  functions.  In  this 
introductory  paper  we  have  used  one-dimensional  examples  f(»r  illustrative 
purposes.  We  will  discuss  applications  in  higher  dimensions  and  extensions 
of  the  theory  of  generalised  distributions  in  following  papers. 
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6  Appendix:  Remainder  Theorem  for  Hy¬ 
perdistributions 

In  this  appendix  we  give  (he  remainder  (heorem  for  h\ perdisirihntions.  I  he 
remainder  formula  is  of  greal  use  in  ralrulating  1  agrange  and  Catirhy  es¬ 
timates  for  the  terms  neglected.  This  feature  of  our  expansion,  like  (he 
presence  of  the  adjustable  ’  radius’'  parameter,  is  unique  to  our  expansion 
and  it  is  not  shared  hy  other  orthogonal  funrti«)ns  expansion. 

The  general  form  of  the  Taylor  expansion  for  two  distinct  points  y  and  z 
is; 

Fiy)^  F{2)  ]  {y-  z)S7F{z)  t  ...  R"*'  (70) 

/?"”  -  0"  V"'’  F{t)/n\  (71) 

The  formula  for  the  remainder  is  in  fact  an  identity  which  is  proven  by 
recursion  integrating  by  parts. 

We  now  give  the  appropriate  remainder  for  the  two  dual  expansions  dis¬ 
cussed  in  section  1. 


6.1  Taylor  (local)  expansion 

We  use 


F{t)  -  x') 

—  6{x'  x) 

y  -  x'  X 

t 

Z  X 


(72) 


We  have 


fx'  X 

6(x'  -  x)  -  ^(x')  xv^(^0  i  *  I  df{x'  X  /)’*  V" ' '  )/”■  (^'^) 

1  hen  we  have,  using  Dirac  ’s  identity 

/(x)  /  ^(x  x')f(x')dx' 

no 


/(O)  xv.f(0)  I  -  ' 


n  I  1 


(74) 


19 


with  the  remarkable  relation 


/  4  <x>  fT*  r 

/  dr' fix')  -  r  /)"  V"'' ^(f)/n! 

J  oc  Jt' 

rcif(T  tr  v"'’/(0/”^ 

Jo 


(75) 


6.2  Moment  (global)  expansion 

We  set 


F{t)  -  i{T  -  x') 

y  -  X  -  x' 

z  ~  X.  (76) 


We  then  have 

S{x  -  x')  -  f{x)  -  t'6{t)  t  ...  -f  j  dt{x  -  x'  -  <)"  V  ^(0/"' 
Using  Dirac  identity  we  conclude 


/(:r) 

/?"'  * 


r”'/(T')6(r  x')dx' 

J  •  oo 

^  4  oc  ^  +00 

I  f{x')dx'  ~  \/S{x)  I  x'f[x')dx'  +  ...4  /?”*’ 

J  -OO  J --00 

dx' fix')  r  ^  di{x  -  x'  -  <)"  V""  ^(0/”! 

J  -  OO  J  T 

/•  +  OO  ^X* 

/  dx'fir')  {-my  x')''v^^'f,ix-y)dy/n\  (78) 

-  c»  JO 


This  remainder  formula  allows  us  to  estimate  correctly  the  errors  when  the 
hyperdistribntion  expansion  is  truncated. 

We  observe  that  the  Rodriguez  expansion  for  a  (Jaussian,  given  in  section 
4  is  also  a  Taylor  expansion  i.e.  (seltingy  -  r  4  A.z  “  x) 


^x(r  )  ^)  ^  v"^i(^)/n! 

fi-n 


2U 


(79) 


Multiplying  both  sides  by  a  good 
obtain 

/+(» 

Sx{t  t  A)<^(x)dr  — 

-  OC 


Taking  the  limit  X  0  of  both  sides  we  conclude  for  the  Taylor  expansion 
of  the  test  function  i.e. 

OO 

^!.{-A)=  ^(-A)"  V"<?i(0)/n!  (81) 

n— 0 

Therefore  the  test  function  miJst  be  real  analytic  to  conclude  that 

6{t  +  A)  -  e^^Six)  f;  A"  V"  /i(x)/n!  (82) 

n^O 

To  show  that  caution  musat  he  exercised  in  treating  hyperdistributions  as 
ordinary  distributions  we  point  out  an  example  of  a  function  which  is  smooth 
but  not  real-analytic. 

<j^{x.)  -  (83) 

The  function  is  infinitely  differentiable  {C°^)  and  tapers  exponentially  at 
infinity.  Nevertheless  it  has  zero  derivatives  at  the  origin  and  it  is  therefore 
not  real-analytic.  As  a  consequence,  the  familiar  Taylor  expansion  of  the  S 
function  is,  strictly  speaking,  a  hyperdistribution. 


(test)  function  (^(x).  and  integrating  ,  we 


^  f  \/''S\{T)<t>{x)dx ln\ 

n-O 

f;(-A)"  f""Sx{T)^”<l>(x)dTln\  (80) 

J  ~00 
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7  Figure  Captions 

•  Fig  1  Schematic  relation  among  the  sections  of  the  paper. 

•  Fig  2  The  Ogura  function  (solid  line)  is  plotter!  together  with  the  5 
terms  (dasher!  line)  and  9  terms  (dot-dashed  line)  Rridrigtiez  expan¬ 
sion  to  the  Ogura  function.  Note  improved  convereence  to  the  Ogura 
function  is  seen  in  the  9  term  expansion.  In  both  expansions,  the  pa¬ 
rameter  A  is  chosen  to  he  1.0. 

•  Fig  3  The  Ogura  function  and  5  term  and  9  term  Rodriguez  expansions 
as  in  Fig  1,  except  with  A  -  0.75.  Note  that  compared  to  the  A  - 
1.0  case,  deviations  between  the  Ogura  functirrn  and  the  Rodriguez 
expansions  are  significantly  larger. 

•  Fig  4  One-sided  moment  plots  of  fi„  ~  /^(T')”/(x')r/T'  for  the  Ogura 
function  (solid  line)  and  its  5  term  Rodtriguez  expansion  (dashed  line) 
with  A  -  1.0.  Note  that  the  's  for  the  Ogura  function  and  the 
Rotriguez  expansion  coverage  as  x  increases  for  n  — 0,2,  and  4  while 

for  Ifi'’  ^  term  Rodriguez  expansion  differs  from  the  Ogura  func¬ 
tion  plot  of  /f6(^) 

•  Fig  5  One-sided  moment  plots,  as  in  Fig  3,  for  the  5  term  Rodriguez 
expansion  of  the  Ogura  ftinction  with  A  —  0.75.  Note  that  the  one-sided 
moments  exhibit  larger  variations  before  converging  to  their  larger  x 
values. 

•  Fig  6  Five  term  Ilermite  expansion  (dashed  line)  to  the  Ogura  function 
(solid  line).  While  convergence  to  the  Ogura  function  is  good  over  the 
domain  in  which  the  Ogura  function  exhibits  significant  values,  outside 
this  domain  rapid  oscillations  (’’whipping  tails”)  are  exhibited. 

•  Fig  7  One-sided  moment  expansions  of  the  5  term  Hermite  approxi¬ 
mation  (dashed  line)  to  the  Ogura  function  (solid  line).  Note  that  all 
of  the  /i„  ’s  computed  for  the  Hermite  expansion  deviate  from  the  //„ 
’s  computed  for  the  Ogura  function  as  x  increases,  indicating  the  fail¬ 
ure  of  the  Hermite  expansion  as  a  global  a|>proximation  to  the  Ogura 
function. 
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•  Fig  8  The  extended  monopole  shown  with  three  valnes  of  the  radius 
parameter  A.  The  graph  shows  that  as  A  -+  0,  we  obtain  a  point  source. 

•  Fig  9  The  extended  dipole 

•  Fig  10  The  extended  qnadrupole 

•  Fig  31  Least  square  error  for  the  Rodriguez  expansion  of  a  gaussian 

•  Fig  12  least  square  error  og  Fig  11  in  expanded  scale 

•  Fig  13  the  rapidly  varying  function  /(r)  =  cosKx  ^  with  K  ^ 
10, A  =  3 

•  Fig  14  Expanded  scale  corresponding  to  Fig  13 

•  Fig  15  The  rapidly  varying  function  with  K  -  4,  A  =  3 

•  Fig  16  The  least  square  error  for  Fig  15 

•  Fig  17  Reconstruction  of  a  T  smeared  with  a  gaussian  filter. 

•  Fig  18  Reconstruction  as  in  Fig  17  with  double  smearing. 

•  Fig  19  Simulation  of  the  CIrossberg-Todororic  neural  network  with  hy¬ 
perdistributions. 
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Figure  1, Relations  among  the  sections  of  this  paper. 
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Hyperdistribution  II:  Two  Dimensional 

Analysis 


We  ronliniie  in  tins  paper  Iho  deve|«>pnien(  of  a  new  technique  for  the  ral- 
rulation  of  ronvolnlion  products  and  llieir  inverses.  Wo  focus  in  this  paper 
on  two-  dimensional  problems.  This  is  accomplished  In  constructing  a  class 
of  singular  “functions’’  of  two  variables  .  hyperdislribulions  of  two  variables, 
that  form  a  closed  algebraic  fiebl  with  the  <  on  volution  prodorl  as  the  mul¬ 
tiplicative  operation.  In  this  paper  we  consider  functions  of  two  variables. 
This  “two  dimensionar'  algebra  can  be  applied  to  image  processing.  Further¬ 
more,  in  this  paper  we  use  the  construction  of  hyperdistribulions  to  obtain  a 
novel  parametric  approximation  method  which  is  the  two-dimensional  ana¬ 
logue  of  the  one-  dimensional  expansion  developed  in  our  previous  paper.  We 
demonstrate  the  use  of  our  approximation  method  with  simple  examples. 
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Introduction 


We  coniimie  to  develop  a  novel  tool  for  the  ralrnlatioTi  of  ronvolution 
products  and  their  in\trses.  We  denote  these  new  mathematical  ohjects  by 
“hyperdistrihutions'’  or  “generalised  ftinctions”  .  These  are  singular  “func¬ 
tions’',  which  we  construct  to  form  an  algebraic  abelian  field  with  the  con¬ 
volution  product  as  the  multiplicative  operation. 

Signal  processing  and  analysis  are  natural  applications  of  hyperdistribu¬ 
tions  in  one  (one  dimensional)  variable.  Image  processing  is  a  corresponding 
application  for  hyperdistributions  in  a  tw’o  dimensional  variable.  Tomogra¬ 
phy  corresponds  to  a  three  dimen.sional  variable  and  the  budding  field  of 
space-time  processing  corresponds  to  four  dimensions. 

The  outline  of  the  present  paper  parallels  the  previous  paper  is  as  follows; 
in  section  1  we  introduce  hyperdistributions  in  two  dimensions  heuristically; 
we  discuss  the  convolution  group  and  derive  a  remainder  theorem.  In  section 

II  we  construct  rigorously  hyperdistributions  by  a  modified  Mermite  poly¬ 
nomial  expansion  in  two  dimensions  and  we  use  the  tools  of  the  Christotfel- 
Darboux  theory  to  obtain  sufficient  conditions  for  convergence.  In  section 

III  we  show  that  Causs’  multipole  expansion  in  two  dimensions  is  obtained 
explicitly  as  a  simple  application  of  the  hyperdistribution  inverse.  Finally,  in 
section  IV  we  expand  a  two  dimensional  gaussian  function  in  terms  of  deriva 
lives  of  a  different  gaussian  to  demonstrate  the  use  of  our  new  parametric 
expansion  and  the  concurrent  minimization  of  error. 

W’e  note  that,  in  effect  rve  achieve  a  method  for  establishing  and  approx¬ 
imating  solutions  of  iniegral  equations  of  ronvolution  form  in  two  variables. 
We  demonstrate  with  examples  that  there  are  cases  for  which  our  method 
is  applicable,  but  Fourier  transform  methods  fail.  Comparing  hyperdistri¬ 
butions  with  Fourier  transforms,  our  method  requires  the  calculation  of  the 
moments  of  the  given  functions  and  of  the  kernel  rather  than  calculation  of 
their  Fourier  coefficients.  This  property  motivates  consideration  of  examples 
for  which  our  method  is  preferable  to  Fourier  transform  techniques.  The 
applications  given  are  motivated  by  image  deconvolution  and  by  the  analysis 
of  turbulent  flow’s. 
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1  Heuristic  definition  of  Hyperdistribution 

We  introduce  a  general  approxiiiiaticni  uliose  integral  properties  arc  the  htcns 
of  interest.  Ovir  apjrroxiniation  displays  in  ronfignralion  spare  the  properties 
of  the  rlas.sica!  moment  generating  expansif)n  for  the  F'oiirier  transform  of 
the  prohahility  distrihnf ion.  Moments,  and  even  shapes,  are  shown  to  he 
captured  well  by  our  expansit»n.  In  additioru  our  expansion  allov^s  ronsider 
ation  of  “functions"  which  are  more  singular  than  temperate  (i.e.  Fourier 
analyzable)  distributions,  but  that  can  be  represented  by  infinite  sequences  of 
distributions.  W'e  call  these  “fiinrtifin.s”  generalised  distributions.  We  show 
that  generalised  distributions  form  the  appropriate  frarneworlr  for  carrying 
out  the  process  of  deconvolution,  in  fact,  we  give  a  straightforward  algo 
rithm  .  the  “TtochnerMartin  algebra",  to  compute  explicitly  the  convolution 
inverse  of  any  generaliserl  distribution.  Applications  of  the  method  to  simple 
optical  deconvolutions  are  shown  to  be  straightforward.  1  hese  applications 
have  been  motivated  originally  by  our  studies  of  optics  and  turbulent  flows. 
We  have  since  become  aware  of  the  much  wider  applicability  of  our  methods. 


1.1  Taylor  and  nioment  expansions 


1  he  theory  of  generalised  distributions  is  built  on  ideas  related  to  Itirac's 
delta  function,  which  is  technically  a  “generalised  function'".  The  Dirac  delta 
function,  MS’),  has  the  properties 


f'ir  ,y)(irdy 


/i(r  )^(  v)d3’dt/ 


f'(x  )d.r  -  1 


(1) 


and 

I  j  y]  (2) 

for  any  suitably  smooth  function,  /(x).  We  term  ('(jualion  (2).  “Dirac  s 
ident'-ty'". 

t  he  Dirac  rlelta  function  is  synuru’l ric  in  its  argument .  i,<'./>(.r  x' .  y  y’) 
^{x'  x.y'  y)  and.  since  th»’  Itirac  della  func  tion  is  a  geiterabsed  func  tion, 
we  may  lavlor  expand  f'ir'  x,y'  y)  about  fr'.j/'). 


3 


This  expression  ariH  similar  rmes  lirlow  are  assnnu’d  1o  hold  under  in1egrali<ni 
Equation  (3)  allows  one  to  ronipiile  a  “loc  al"  approximation  to  /(r.y),  since 
if  we  suhstifiite  this  ex|>ansioii  into  Dirac's  identity,  we  recover  the  usual 
Paylor  expansion  of  fir,y)  about  (r.j/)  (O.D). 

f(T,y}  /(O.O)  .  rVJ(n,0)  t  yVJ(n,n)  I  ...  1  (1) 

This  approximation  is  local  in  the  sense  that  it  requires  derivatives  of  /(t) 
at  a  single  point  .r  -  0.  and  in  general  has  a  limited  radius  of  convergence. 
On  the  other  hand,  if  we  expand  /'(yr  T',y  y')  about  r,  we  have 

t',v  y*)  ~  ■  x'V,/i(T,y)  <  ...  (  //”*’  (5) 

When  this  series  is  substituted  into  Dirac’s  identity,  we  obtain 

/(x,y)  Arfi(T)S(y)-  Af'VAlT)^(y)  >  ...  I  («) 

The  coefTicients  A/",  defined  by 


A/"  - 

f  f 

All  - 

1  1  Tf{T,y)<iT(iy 

K  - 

j  j  yf{^-,y)A^<iy 

(7) 

are  simply  the  moments  of  the  fiinefion  f{T,y).  1  herefore.  equation  (6)  is 
an  approximation  of  f{T,y}  involving  global  information  about  /(x,y),  that 
is.  the  moments  of  the  function. 

This  then  may  be  taken  as  a  motivation  for  our  definition  of  a  generalised 
distribution  as  a  function  which  may  be  written  in  the  form. 

OO 

/(x.y)  - 

«-n 

no  '  On 

ni  0  V  -  n|:r^T  » 

—  n2TxV‘  1  2n2Ty^  ,,  *  ”2vy^y  (^) 

The  n,|  values  are  coefficients  given  by 

o,.  (  irA/'-./n!  (fl) 

Note  that  equation  (8)  is  equivalent  to  the  familiar  moment  generating  ex 
pansion  of  probability  theory  (see  equation  b). 
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1.2  The  Convolution  Group 

If  we  have  two  goneraliscfl  <V\st rihntions  /j  and  fi.  Oteir  linear  rornbinalioii 
■^/i  1  f‘f2  siso  a  generalised  distrihiHion  (where  A  and  //  are  real  eoefH 
eients).  The  pth  derivative  of  a  generalised  dislrilnition,V'’/(j', y),  is  a  gen 
eralised  distriliiition.  Also  .  the  ronvolntion  of  two  generalised  distribution 
/i  *  /s  is  a  generalised  distribution.  These  may  all  be  thought  of  as  “closure 
properties  ’  of  generalised  dist ributi<*ns. 

CJeneralised  dislributif)ns  allow  tjs  to  make  effeetive  computation  of  the 
convolution  inverse,  (liven  a  generalised  distribution  /,  the  desired  convolu¬ 
tion  inverse  /nj/j  satifies 


/  * /n|/l  -  «{3-,y)  (10) 

Here  S{T,y)  represents  Dirac's  delta  function,  which  is  the  identity  of  the 
convolution  operation.  VVe  shall  show  by  our  construction  that  In\f\  is  a 
generalised  distribution.  Writing  the  convolution  explicitly, 

f*fn\f\-  I  I  x'.y  ■  y)  r-  f>{T)h{y)  (11) 

It  is  now  necessary  to  compute  the  product  of  the  sums  and  match  coef¬ 
ficients.  Taking  /(x.y)  -  n  «n  G  V"^(3’)<S(y),and  y(r,y)  - 
V"f(x)^(y),  we  see  that  the  computation  of  the  convolution  inverse  is  effec¬ 
tively  the  determination  of  a  collection  of  h„  values,  given  a  set  of  a„  vahies. 
Substituting  into  equation  (11), 

V"^(3-.y)  *  V’"^(r,y)  -  V’-’"-^(T,y)  (12) 

This  gives 

r  9 

or  equivalently  with  r  p  I  (/ 

no  r» 

1.]  X! ,,V’’('i{T.y)  A(.r.y)  (11) 

P  n  r  n 

Matching  coefTicients  on  the  left  and  right  hat\d  si(h-s  implies  that  only  tin' 
(1  term  survives.  The  result  is  a  linear  system  <»f  equations  for  the  h 


r 


values  in  terms  of  the  a  values.  It  is  easier  to  see  the  behavior  by  'vrifiiip;  the 
first  equations  in  this  linear  system. 


00^*00  I 

Qoo^in  ‘  '^10^11(1  b 

O0()^>0l  *  <^01^00  ~  b  (1'^) 

and  so  forth.  Thus,  we  ran  see  the  computation  of  the  terms  of  the  ronvolu 
tion  inver.se. 


1.3  Fourier  Transforms  and  Simple  Examples 

Lastly,  we  shall  want  to  consider  the  Fourier  transform  of  a  generalised  dis 
tribution.  Again  taking  /(j’,y)  -  o  ?/)'  immmecliately 

evaluate  the  Fourier  transform  as; 

I  ik)"  (Ifi) 

J  n--0 

where  for  the  last  step  we  have  used  ■  (lA’)''.  fhe  Fourier 

transform  of  a  generalised  distribution  may  seen  to  be  a  power  series. 

The  Fourier  transform  of  a  function  is  the  “moment  generating  function" 
because  of  equation  (9).  Thus  the  basic  requirement  for  the  validity  of  writing 
a  function  as  a  generalised  distribution  is  that  tis  moments  he  finite,  or  more 
stringently,  that  its  Fourier  transform  be  real  analytic. 

Simple  explicit  examples  of  convolution  inverses  are  obtained  from  the 
standard  (lieen's  functions.  The  two  dimensional  Helmholtz  equation,  for 
example,  can  be  treated  equally  easily  with  generalised  distributions  and  with 
Fourier  transforms  to  obtain  both  the  Green’s  function  and  its  convolution 
inverse.  Thus  ,  we  have 

(V^  i  )/1)r  (17) 

Ap])lication  of  equation  (I'i)  slu>w’s  that  only  and  hi  contribute  to  (he 
convolution  inverse 

hi\<  ?  I  fi{r.y)/2  b  (x,y) 


(18) 


•••  • 


1  liis  rrsnlt  is  f'asily  rrprndiirrf)  l)y  using  Funrirr  Irnnsfornis  anti  is  easily 
verified  using  equation  (10). 

In  eunfrasi,  flie  eonvolntion  inverse  r>f  the  (lanssian  fiinrtioji  cannot  he 
obtained  as  a  Fourier  Iransfonri.  but  it  is  easily  obtained  using  generalised 
distributions.  It  may  be  verified  that  for  a  (lanssian  fiinelif)n  of  width  n  ,  the 
convolution  inverse  is  readily  obtained  from  equation  (lb),  while  calculation 
using  Fourier  transforms  leads  to  a  divergent  result,  ('alrulation  of  mornnets 
and  use  equation  (15)  yields 

<I2t.  '  (fT/\/2)"/rj! 

-  (  l)"(e/V2)"/»!  (19) 
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2  Rigorous  Construction  of  Hyperdistribu¬ 
tions 

In  Ibis  sorlidH  we  first  modify  tfic  classical  expansion  in  the  llermife  polyno 
rnials.  This  expansion  is  “superior”  lo  a  power  series  expansion  in  that  the 
terms  are  orthogonal,  making  the  error  orthogonal  to  the  approximation.  Onr 
modification  maintains  this  advantage.  Fnrthnremore,  the  standard  Christof 
fel  -Darhonx  analysis  gives  a  very  useful  sufficient  condition  for  convergence 
that  onr  expansion  inherits  from  the  Ilermite  polynomial  expansion.  Since 
onr  expansion  utilises  the  Rodriguez  formula,  for  the  Hermite  polynomials, 
we  call  it  the  “  Rodriguez  expansion’^ 


2.1  Modified  Hermite  Expansion:  Rodriguez  expan¬ 
sion 

In  effect,  we  discuss  here  a  systematic  pointwise  approximation  of  generalised 
distributions.  This  approximation  method  is  analogous  to  the  method  de¬ 
veloped  by  Temple  lo  approximate  distributions  by  smooth  functions.  The 
Temple  method  has  become  widely  known  through  Lighthill's  monograph: 
Fourier  Series  and  Generalised  Functions. 

Generalised  distributions  are  approximate  by  a  modification  of  the  classi¬ 
cal  expansion  in  flermite  polynomials  (hereafter  called  the  “Hermite  expan¬ 
sion’’)  which  is  suggested  by  the  Rodriguez  formula; 

(  l)"'"‘//„(x)fC(T/)r  (20) 


To  see  how  the  Hermite  expansion  can  be  transformed  into  an  approximation 
for  generalised  distributions,  consider  a  function  /(x,y)  which  we  want  to 
represent  as  a  generalised  distribution.  Multiply  f{x,y)  by 


*■  .tr 

r 

ttA 


(2!) 


then  expand  the  resulting  expression  in  terms  of  the  “scaled  Hermite  poly 
nomials” 


//,^(x)  A  "//..(x/A) 


(22) 
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w}k)Sc  (Icfiriil ion  is  jnstifird  l)y  ihr  fronnila  (2r>)  hrlow.  \\c  fihlain  llie  rx 
pression 


MtiUiply  both  sides  of  the  equation  by  the  normalised  Caiisstan  introduced 
by  Temple,  i.e. 

hix.y)-"'  'I  (2-1) 

TT  A 

And,  observe  that  by  rescaling  the  Rodriguez  formula  (20)  we  can  write 


f(x,y) 


Li: 

n"  0  m  "  ( 


(-  l)"'-W„^(x)//^(t/)Mx,y)  -  ^:V;'6,{r,y)  (25) 

The  final  result  is 

/(»■,!<)  -  £  2;  (20 

n  -  0  tn  -  0 

where  we  have  introduced  the  coefficients 


-  (  1) 


n4  m  . 


n,m 

\  2ri  t  2*n  - 


a 


n,m 


When  w'e  let  A  fend  to  zero,  we  have  a  representation  of  /{r)  as  a  generalised 
distribution.  We  call  the  expansion  (26)  the  “Rodriguez  expansion  for  /(x) 
with  width  A”.  The  width  parameter  is  a  novel  feature  of  the  Rodriguez 
expansion  when  compared  to  standard  expansions  in  complete  sets  of  basis 
functions.  The  standard  expansions  do  not  contain  a  free  parameter.  The 
advantages  of  the  Rodriguez  expan.sion  will  be  demonstrated  below  in  the 
context  of  the  theory  of  generalised  distributions  and  of  convolution  inverses. 

We  consider  as  an  example  a  function  familiar  from  the  analysis  of  tiir 
bulenre  spectra,  the  “Ogiira”  function 

/(t,.v)  ■  Ae  (28) 


With  A'  determined  by  the  normalisation  condition, 

I  I  f{x.y)dTfiy  I 


(29) 


9 


we  have 


TV  =  (l/4)r*(l/4)  =  (1.8128)*  (30) 

The  Hermite  expansion  exhibits  pointwise  convergence,  albeit  with  “whip¬ 
ping  tails”.  By  contrast,  the  Rodriguez  expansion  exhibits  pointwise  conver¬ 
gence  without  “whipping  tails”.  As  a  consequence,  the  Rodriguez  expansion 
also  yields  a  global  approximation  by  accurately  representing  the  moments. 
Furthurmore,  the  Rodriguez  approximations  has  an  adjustable  parameter 
which  may  be  selected  to  optimise  the  rate  of  convergence.  The  Christoffel- 
Darboux  theory  gives  the  necessary  condition  for  convergence  in  appropriate 
£*.  For  the  Hermite  expansion,  we  have 

jj  e~^^'^'^'^p{x,y)dxdy  <  oo  (31) 

Inserting  the  transformation  that  leads  to  the  Rodriguez  expansion,  we  find 

j  j  e**^‘^/*(®,y)dxdy  <  oo  (32) 

for  the  corresponding  £*  definition  for  the  Rodriguez  series. 

2.2  Definition  of  Hyperdistributions 

In  order  to  describe  the  process  of  antidiffusion,  we  have  introduced  a  class 
of  highly  singular  “functions”,  that  are  precisely  defined  in  this  section.  Our 
process  for  defining  hyperdistributions  parallels  the  Temple  defiii’fion  (gen¬ 
eralised  function)  as  a  good  sequence  of  good  functions.  Good  functions  are 
smooth  and  tapered.  More  precisely,  they  are  total  point  functions  that  a’-e 
differentiable  to  all  orders  (0°°),  and  decay  at  ±oo  faster  that  any  power. 
Good  functions  play  the  role  of  “testing”  a  sequence  of  good  functions  for 
weak  convergence.  In  fact,  a  sequence  of  good  functions  is  a  distribution  if 

lira  /  /  (f>{x,y)f„{x,y)dxdy  <  cx>  (33) 

n—oo  J  y_o(, 

for  all  good  functions  4>. 

Since  we  conceive  of  hyperdistributions  as  “generalised”  distributions,  we 
are  in  fact  implementing  a  second  order  generalisation  of  functions.  Con¬ 
sequently,  we  need  a  double  test  as  a  convergence  criterion.  We  implement 
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this  criterion  by  introducing  very  jood  functions  G\{x,y)  with  the  following 
properties: 

1.  G\{x,y)  is  smooth,  that  is,  differentiable  to  all  orders,  C°°. 

2.  G/t,{x,y)  is  essentially  compact,  i.e.  it  has  a  gaussian  decay  at  ±oo  : 

We  will  assume  for  convenience  that  G\  is  normalised  to  unity: 

j  j  GK{x,y)dxdy  =  1  (34) 

We  define  the  width  of  Ga  by 

A*/4  =  /  {x  -  x)'^{y  -  yYGA{x,y)dxdy  (35) 

J  —  OO 

A  primary  example  of  very  good  function  is  the  gaussian,  which  we  denote 

by  h{x,y)'. 

=  ^--^2  (36) 

We  now  introduce  a  sequence  of  very  good  functions  defined  by 

W\Ut)  =  l:H||®VllM5,t)  (37) 

11=' 


The  sequence  where  lambda  is  a  nonnegative  real  and  n  is  a  natural 

number,  is  a  good  sequence  if,  for  all  good  functions  (f>  and  for  all  very  good 
functions  Ga,  there  exits  a  Ao  such  that,  for  all  A  >  Ao, 


lim 

n— *Qo 

A-.0 


<  OO 


(38) 


We  note  that  e*‘^*^(x,j/)  are  hyperdistributions.  The  sum  (hyperdistribu¬ 
tion) 

f^au<SVH{x,y)  (39) 

*=o 
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mil  thus  1)0  virnvrd  cifhrr  as  a  srijurTicr  nf  “pood"  disi rihul ions  as  A  >  (I: 


G' (10) 

*■  (» 

or,  as  7?  >  oc.  as  sriiurnre  of  good  fiinrtioiis: 

OC 

(4!) 

Ir-n 

The  latter  representation  is  a  Rniirifjufz  rrpnvsion.  The  Rodriguez  formula 
for  the  Hermite  polynomials  can  he  used  to  show  the  derivatives  of  a  gaussian 
form  a  complete  set  of  orthogonal  polynomials  in  an  spare.  And  thus  the 
Rodriguez  expansion  yields  a  very  useful  point  function  approximation  to 
any  hyperdistrihution: 

OO 

X  »»,( *  1  )‘"'Ht(x/\)IMy/X)S,{x,y)IX’' ' '  (U) 

k- 0 

where  //„(.r)  denotes  the  Hermite  Polynomial  in  t  of  order  n. 
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3  The  Multipole  Expansion 

\Vf  slarl  wilti  flic  f;ni)ilinr  I’oissoti  rqnniioii  <ff  j)t)1t‘niial  tlirnry, 


VV  P  (13) 

Tliis  equation  is  rewritten  with  the  lu'lp  of  the  (infinite  domain)  (ireen's 
function 


(!{T,y)  /u(?)/2/Tr 

V'a  f{x,y) 


o  7 

r  -  3-  1  y‘ 

(11) 

We  can  thcfi  rewrite  the  “potential”,  <^.  in  terms  of  the 
tion’',  p  as 

“charged  dislrihu 

4)  -  G  *  p 

(15) 

Introduce  Q  with  the  property 

O  *  p  fi 

(If,) 

Convolving  hoth  sides  of  eq(‘15)  with  Q  and  using  the  rommutalive  and 
associative  properties  of  the  *  product,  we  find 


Q  *  if)  —  Q  *  ((7  *  p) 

=  {Q  *  p)*  C 

-  (?  (IT) 

Solve  eq(13)  for  (f>  in  terms  of  the  given  G  by  computing  the  convolution 
inverse, 

OO 

(48) 

fro 

which  is  Causs'  rntiltifjole  series  with  coefficients 

(  ^^) 

where  (g)  denotes  tensor  product  and  r**"  is  the  tensor  power  of  the  3  vector 

X. 


13 


Suhsfil  ntiiig  (  1!))  into  (  IS)  we  find  ;i  fnttiilinr  rxprrssifiit 

'/’(■»'•?/)  X]  ^  ^  i  f  pi  r,y)d.yHy)  ^  (5(1) 

k  n  '•■•  ■’  ■'  '  27r 

wliicli  is  a  statuiard  rrsnit  in  potrntial  (linorv. 

\\r  ran  intrrpr('1  ()ic  Rodripiic?  rxpansioti  as  a  gnu'ralisrd  niiiltipole  rx- 
pansion  Ilia)  intlndcs  a  si7<'  or  ratlins  paranilrr.  A.  rhiis  oiir  inonopolr 
generatps  the  point  sonrre.  wliirh  is  a  Dirar  drlta-fntutit)?!  as  a  ganssian  of 
width  A.  We  recover  the  standard  innltipoles  wlien  A  >  0  .  Our  expansir>ii 
has  the  form 


/(r.i/)  -  f  ^  I  ...  (51) 

The  ingredients  of  the  expansion  are  the  basis  functions 

and  tlie  coefricients 

all  of  which  depend  on  (he  size  parameter  A.  The  basis  functions  are,  explic¬ 
itly 


^A(T,y) 


e 


7rA2 


-  2t 


'^y^A(x,y)  -  -2y 


e  ^ 

----- 

e  ^ 


The  corresponding  expansion  coefficients  are  given  by 


«o  ~  (A  >  0) 


(51) 

(55) 

(56) 

(57) 


This  coefficient  represents  the  total  area,  or  alternatively,  the  total  mass  or 
charge  of  the  source  function.  Also 


f  f  ^  ^ 

I  I  ^  ^f{3-,yH3’dy,  (A  >  0) 


(58) 
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witicli  rr[)TTSf'iits  t  hr  iH'l  riipolr  of  f  hr  soiirrr.  In  iiirchafiiral  Irnns  drfinrr 
Ihr  rrnlrr  of  mass  of  thr  sourer  fiinrlion.  I  inatly 

~  /  I  ^'f{3-,y)drtiy.  (A  >0)  (59) 

which  represents  the  qnadnipole  of  the  source  distrihiitions.  In  mechanical 
terms  this  is  the  momrni  of  inertia  of  the  mass  distribution.  Higher  order 
basis  function  and  coefficients  have  analogous  interpretations. 


4  Expansion  of  a  Function  in  Rodriguez  Se¬ 
ries  and  its  Optimization 

Our  expansion  preserves  tiie  dassieal  |)rf)perlies  (liat  are  (lerival)le  for  or 
ihogonal  polynomial  expansions  (as  opposed  to  power  series),  hut  also  adds 
an  important  new  feature;  the  si7e  or  radius  parameter  that  generates  the 
CJaussian  pieliire  of  point  niono|>ole.  <li})oie.  piiadruprde.  etr,.  to  a  seeiiario 
in  which  we  can  allow  ffir  ext(Mide<l  sources.  W'e  find  frotrt  the  convergcTice 
condition  dt-rived  from  the  ChristofTel  Darhoux  theory  in  sectif)n  II  that  con¬ 
vergence  holds  generally  on  a  semi-infinite  range  of  A.  'I  his  remarkairle  free 
dom  is  ex|)loited  in  this  section  tf»  optimiz<‘  Mie  rate  of  convergciicf'  of  tfie 
expansion.  I  his  optimization  results  in  determining  the  value  of  A  for  which 
a  i  imum  numher  of  terms  is  determined  in  order  to  obtain  a  given  toler 
arue.  We  measure  the  tolerance,  as  usual  .  hy  the  least  square  fit  integrated 
over  the  entire  function.  W'e  define  the  error 

r(.V.A)  //  dTdy{f{x,y)  (fiO) 

•'  •'  k  0 

W'e  apply  the  formula  to  a  standard  gaussian,  i.e.  we  take 

(61) 


we  then  find,  after  some  calculations 


(  If 2*^  2’‘(2r7  -  2A-  1)!! 

k'!(2n  -  2A’)! 


t  (  1)''2  ^"/n! 


where  the  double  factorial  contains  rrrdy  od<l  terms.  W’e  can  now  express  tlie 
standard  gaussian  in  the  Ilodriguez  form; 

(a-’4d)  =»  /.  -210 

A  profif  based  on  useful  linear  operator  relations  is  as  follows.  We  ran  check 
bv  differentiation  in  f  and  .r  that 


H^-y) 


\\  r  nf>\v  Id 

As  ;i  cotisftjiirnci'.  wt-  cnii  wrilc 

7r  A- 

\\{'  !i(nv  coiisidcr  fin'  idrnlitv 


1/  A' 


*■  r' , 


expatiHitig  (Itr  It.ll.S. 


(  »• 


ttA’ 


(1  A^)"  ,, 

2^  2?n„!  ~ 


ti  0 


I  hr  !nr1h«Kl  of  proof  will  lie  ttsrd  lalcr.  Ftirt hnrtiioro,  it  prrnirirs 
cfiock  on  t hr  ratlirr  (lifTirnlt  calrnlations  of  corfnrirnls,  \Vr  sre  hy 
that  A  I  is  optimal.  A^’r  also  ronsirlrr  (hr  function 


/( T,y)  -  ro.s( kx )ros{ky ) 


ttA^ 


fti",) 


(fir.) 


(fi?) 


(68) 

a  welcome 
insprcl  ion 


(69) 
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5  Conclusions 

\\  (■  }tavc  scrii  lhal  E<'iH'rnlis<‘il  rilml  ion  pr<>\i(lc  a  iiu'fliixi  for  ‘.olvinn 
I'rt'dliohii  inlcgra!  ('(juatimis  of  tlio  rmi volii) inn  1\i>i'  wtiicli  is  <i>tnj)cf itivr 
witli  Uk'  rrirrenl  rnrllmds  iliat  niiploy  Fourier  transfonns.  V\'c  have  aiso  seen 
Oiat  generalised  fiisfril)ii1if>ns  ran  lie  approximated  nnnierieallv  1»\  sctpienees 
of  srnootl)  functions.  l  ids  pr<icerli)r«-  is  analogous  to  titat  of  a|)|)roxintatiTig 
Dirac  tirlin  fiinctif>ns  liv  sequences  f>f  narrowing  (laiissian  functions.  In  tins 
secoiul  pafier  we  liaA-e  used  t \vf)  dimensional  examples  for  illnst  rati  ve  purposes 
W  e  will  disc  uss  applic  ations  in  higher  dimensions  and  extcmsions  of  t  lie 
theory  cif  peneralised  disf rihul ic>ns  separafcdy. 
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6  Appendix:  Remainder  Theorem  for  Hy¬ 
perdistributions 

hi  this  npix'iulix  wr  f'ivr  thr  rfinaiiuh'r  llirori'in  for  liypf’rdisl  ri  but  ions.  I  iic 
mnniiidrr  forinn1;i  is  of  ^rral  use  in  ralciilntinR  I.agrangr  anri  (’ancliy  fs 
tiinate's  for  tlic  ti'rins  nr^lrctnil.  I  his  h’atnrn  of  onr  expansion,  like  (he 
presenre  of  (he  adjiisiahlc  ''radins"  paranH-ler.  is  iiniipn’  to  onr  expansion 
and  it  is  not  shared  Ity  other  orthogonal  fnneiions  expansion. 

I  he  general  form  of  t  he 'i  avion  expansion  for  t  wo  dlst inct  points  j/ and  r 
is; 

Fiy)  Viz-)  \  [y  z)  ■  V  l'(z^  ^  >  IF  ‘  '  (7ll) 

/)"  '/'{/■)/»'  (71) 

1  lie  formula  ff)r  the  remaind"r  is  in  fart  an  identity  which  is  proven  by 
renirsion  integration  by  parts. 

We  now  give  the  afipropriate  remainder  for  the  two  dna!  expansions  dis 
enssed  in  section  I. 


6.1  Taylor  (local)  expansion 

We  use 


f  (r)  f'{T  t') 

y  ^ 

z  t' 


We  have 
fi{x'  t)  - 


riien  we  liave,  using  Dirac's  identity 


(72) 


/)”  «  (73) 


/(.f)  /  /'(.r  T')/(;r')f/.r' 

./  li” 

/(tt)  .  v/(n)  I  ...  I  /r" 


(H) 


t 


wilti  l!\r  mnarkahir  relation 


r”  '/(/')  n' 


(7-.) 


6.2  Moment  (global)  expansion 

Wo  set 


/I,.)  >.1.:  y>\ 

V 


(Tfi) 


We  tlieii  liave 


f{T  r')  o'  *  ... 

I/siiig  Dirae  identity  we  conclude 


V’7(/);r,!  (77) 


/(•'•)  Hr)  I  !(T')d"o'  (78) 

'fin's  remainder  foriinila  allows  ns  to  estimate  correctly  tlie  errors  when  the 
hy{)erdist rihnt ion  expansirjn  is  truncated. 

We  observe  that  the  Hodrignez  expansion  for  a  (Janssian.  given  in  section 
1  is  also  a  faylor  expansion  i.e.  (setting  y  ~  or  !  A,  z  -  x) 


•  A)  ^  V'’^,(T)/n!  (79) 

n  -  0 


Multiplying  both  sides  by  a  good  (test)  function  <f>{x),  and  integrating,  we 
obt  ain 

/  ^A(.f  i  ii)</.(f)dr  -  y(  A)"  5?/  fHr)V’'d{T)dr/n\  (80) 
J  ./  /?" 

'faking  the  limit  A  »  0  of  botli  sides  we  conclude  for  the  faylor  expansion 
of  the  test  function  i.e. 

d{  V)  y{  A)"(??VX0)/n!  (81) 

T»  0 
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(  fnT(T'ir('  till-  1r<t  iinist  fii-  H  ill  iitiith  lii  In  (oiulildi’  tliiit 

X 

f'ir  ■  <,)  V(_\r r”^fr)  »'  fv:>) 

•  r  fl 

'{(V  slifiw  (lint  rntitinti  tTiiict  lu-  rxcrcisf'd  in  lr<’;t(iiie  inin'rdi'^t rilint inns  ns 
oritinnrv  (iislriliiitions  wf  point  out  nn  rx.'irnplo  of  n  fnni  tioii  Indt  is  smooili 
hill  not  rcni  nnnlyru  , 

dir.i,)  ,  '  (q;t) 

The  function  (t>  is  infinitely  {lifferent inlile  (f'"^  }  nnrl  tnpers  exponent inlly  at 
infinity,  Nf'vert lieless.  it  lias  zero  rierivatives  at  the  nrij^in  and  it  is  therefore 
not  real  analyt ic.  As  a  consequence,  the  fainilar  daylor  exjiansion  of  the  ^ 
function  is,  .strictly  speaking,  a  livpcrdistrihntion. 


